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Abstract 

A third-order Energy Stable Weighted Essentially Aon- Oscillatory (ESWENO) fi- 
nite difference scheme developed by Yamaleev and Carpenter [1] was proven to 
be stable in the energy norm for both continuous and discontinuous solutions of 
systems of linear hyperbolic equations. Herein, a systematic approach is presented 
that enables “energy stable” modifications for existing WENO schemes of any order. 
The technique is demonstrated by developing a one-parameter family of fifth-order 
upwind-biased ESWENO schemes; ESWENO schemes up to eighth order are pre- 
sented in the appendix. New weight functions are also developed that provide (1) 
formal consistency, (2) much faster convergence for smooth solutions with an ar- 
bitrary number of vanishing derivatives, and (3) improved resolution near strong 
discontinuities. 

Key words: high-order finite difference methods, weighted essentially 
nonoscillatory schemes, energy estimate, numerical stability, artificial dissipation. 


1 Introduction 


A new, third-order weighted essentially nonoscillatory scheme (called Energy- 
Stable WENO) has recently been proposed and developed by the authors of 
the present paper [1], In reference [1], we prove that the third-order ESWENO 
scheme is energy stable, that is, stable in an L 2 energy norm, for systems 
of linear hyperbolic equations with both continuous and discontinuous solu- 
tions. Stability is explicitly achieved (by construction) by requiring that the 
ESWENO scheme satisfies a nonlinear summation- by-parts (SBP) condition 
at each instant in time. Thus, L 2 strict stability is attained without the need 
for a total variation bounded (TVB) flux reconstruction or a large-time-step 





constraint [2], [3] and [4], Herein, we generalize and extend the third-order 
ESWENO methodology [1] to an arbitrary order of accuracy. Similar to the 
third-order ESWENO scheme, the new families of higher order (up to eighth 
order) ESWENO schemes are provably stable in the energy norm and retain 
the underlying WENO characteristics of the background schemes. Numerical 
experiments demonstrate that the new family of ESWENO schemes provides 
the design order of accuracy for smooth problems and delivers stable essen- 
tially nonoscillatory solutions for problems with strong discontinuities. 

Another issue that is addressed in this paper is the consistency of the new 
class of ESWENO schemes. The consistency of any WENO-type scheme fully 
depends on a proper choice of the weight functions. On one hand, for smooth 
solutions the weights should provide a rapid convergence of the WENO scheme 
to the corresponding underlying linear scheme. On the other hand, the weights 
should effectively bias the stencil away from strong discontinuities. The high- 
order upwind-biased WENO schemes with conventional smoothness indicators 
that are presented in reference [2] are too dissipative for solving problems with 
a large amount of structure in the smooth part of the solution, such as direct 
numerical simulations of turbulence, or aeroacoustics[5], [6]. Furthermore, as 
has been shown in references [7] and [8], the classical weight functions of the 
fifth-order WENO scheme fail to provide the design order of convergence near 
smooth extrema, where the first derivative of the solution becomes equal to 
zero. New approaches are proposed in references [7] and [8] to improve the 
error convergence near the critical points. Although these new weight func- 
tions recover the fifth order of convergence of the WENO scheme near smooth 
extrema, the problem persists if the Erst- and second-order derivatives vanish 
simultaneously [8]. An attempt to resolve this loss of accuracy is presented in 
reference [8]. This proposed resolution provides only a partial remedy for the 
problem; the same degeneration in the order of convergence occurs if at least 
the Erst three derivatives become equal to zero. To fully resolve this problem, 
we propose new weights to provide faster error convergence than those pre- 
sented in reference [8], and impose some constraints on the weight parameters 
to guarantee that the WENO and ESWENO schemes are design-order accu- 
rate for sufficiently smooth solutions with an arbitrary number of vanishing 
derivatives. 

This paper is organized as follows. In section 2, energy estimates for the con- 
tinuous and corresponding discrete wave equations are presented. In section 
3, we present a one-parameter family of fifth-order WENO schemes; one value 
of the parameter yields a central scheme that converges with sixth-order ac- 
curacy. In section 4, we present a systematic methodology for constructing 
ESWENO schemes of any order and demonstrate the methodology by trans- 
forming the family of WENO schemes presented in section 3, into a family 
of fifth-order ESWENO schemes. In section 5, we analyze the consistency of 
the new class of ESWENO schemes and we derive sufficient conditions for the 
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weights functions that ensure that the ESWENO schemes are design-order ac- 
curate regardless of the number of vanishing derivatives in the solution. The 
tuning parameters in the weight functions are also optimized. In section 6, we 
present numerical experiments that corroborate our theoretical results. We 
summarize and draw conclusions in section 7. 


2 Energy Estimates 


Consider a linear, scalar wave equation 

f + f = 0, f = au, t> 0, 0 < x < 1, 

u( 0, x) = Uq(x) 


where a is a constant and u 0 (x) is a bounded piecewise continuous function. 
Without loss of generality, assume that a > 0, and further assume that the 
problem is periodic on the interval 0 < x < 1. Applying the energy method 
to equation (1) leads to 


d_ 

dt 



0 


( 2 ) 


where || ■ \\l 2 is the continuous L 2 norm. Thus, the continuous problem defined 
in equation (1) is neutrally stable. 

We now develop using mimetic techniques (see [1] or [9]) a class of discrete 
spatial operators that is neutrally stable or dissipative. The continuous target 
operator used for this development is the following singular perturbed wave 
equation: 


N 

f +S= <>o. o<*<i, 

n = 1 x 7 

u{ 0,x) = Uq(x) , 


( 3 ) 


where // = n{u) is a non-negative C°° function of u. As before, we assume 
that equation (3) is subject to periodic boundary conditions. Our goal is to 
match each spatial term in equation (3) with an equivalent discrete term that 
maintains neutral stability (or dissipates) of the discrete energy norm. 

We begin by showing that the terms on the right-hand side of equation (3) 
are dissipative thereby ensuring stability. Multiplying equation (3) by u and 
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integrating it over the entire domain yields 


Id 2 
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-L O 
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i 

o 



( 4 ) 


Integrating each term on the right-hand side by parts and accounting for 
periodic boundary conditions yields the following energy estimate: 


d_ 

dt 




( 5 ) 


All the perturbation terms included in equation (3) provide dissipation of 
energy. 


Turning now to the discrete case, we define a uniform grid Xj = j Ax, j = 0, J, 
with Ax = 1/J. On this grid, we define a flux f = au and its deriva- 
tive f ;c = au x , where u = [u(x 0 ,t), . . . ,u(xj,t)] T and = [u x (x 0 , t), , 

u x (xj,t)] T are projections of the continuous solution and its derivative onto 
the computational grid. Next, we define a pth-order approximation for the 
first-order derivative term in equation (1) as 

or 

— = Df + 0(Ax p ) . (6) 


Placing a mild restriction on the generality of the derivative operator (see [1] 
or [9]), the matrix D can be expressed in the following form: 

D = P~ 1 [Q + R\ ; Q + Q T = 0 

R = R T ; v T i?v > 0 (7) 

P = P T ; v T Pv > 0 

for any real vector v ^ 0. By choosing the matrix R as a discrete analog of 
the dissipation operator in equation (3), we have 

N 

R=Y. TVAprf , (8) 

n = 1 
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where A is a diagonal positive semidefinite matrix and D\ is the difference 
matrix 


0 


D 1 = 


-1 1 


0 


By using the SBP operators [eqs. (6)-(8)], the semi-discrete counterpart of 
equation (1) becomes 


<9u 

~dt 


+ P~ 1 Qi 


N 

Y P _1 Di n A[Di n ] T f j 

n= 1 


( 9 ) 


where f = au, u = [u 0 (t), Ui(t ), . . . , Uj(t), ] T is the discrete approximation 
of the solution u of equation (1), and Q and A are nonlinear matrices (i.e., 
Q = Q(u) and A = A(u). To show that the above f ini te- dif ference scheme is 
stable, the energy method is used. Multiplying equation (9) with u T P yields 

i j v „ 

--|| u ||^ + a u ^ u =-aE([^r] Tu ) - / V[T> 1 "] t u, (10) 


where || • ||p is the P norm (i.e., ||u||p = u r Pu). Adding equation (10) to its 
transpose yields 


d_ 

dt 


u||p + au 1 


(Q + Q T ) u = 


N T 

-2aY([Di n ] T u) AlD^fu. 

n= 1 


in) 


If we account for periodic boundary conditions and the skew-symmetry of Q, 
then the second term on the left-hand side vanishes, and the energy estimate 
becomes 

d N rp 

-||u|| 2 P = -2«£([Br] T d A[ZVfu<0. (12) 

al n= 1 


The right-hand side of equation (12) is nonpositive because the diagonal ma- 
trix A is positive semidefinite [v T Av > 0 for all real v of length (J + 1)] and 
a > 0; thus the stability of the finite-difference scheme given by equation (9) 
is assured. This result can be summarized in the following theorem: 

Theorem 1 The approximation [eq. (9)] of the problem [eq. (1)] is stable if 
equations [(6)-(8)] hold. 

Remark 2 Despite the fact that the initial boundary value problem [eq. (1)] 
is linear, the finite-difference scheme [eq. (9)] constructed for approximation 
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of equation (1) is nonlinear, because the matrices Q and A (and in principle 
P) are assumed to depend on the discrete solution u. 

Remark 3 The only constraints that are imposed on the matrix Q and the 
diagonal matrix A are the skew-symmetry of the former and positive semidef- 
initeness of the later. No other assumptions have been made about a specific 
form of the matrices Q and A to guarantee the stability of the finite-difference 
scheme [eq. (9)]. 

Remark 4 The discrete operators that are defined by equations [(6)— (8)] are 
similar in form to those that are used for conventional SBP operators (See refs. 
[9,10]). What is new, however, is the fact that the matrices Q and A depend 
on u. 
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Fig. 1. Extended six-point stencil S' 6 , and corresponding candidate stencils Sll , Sl, 
Sr, and Srr for one-parameter family of fifth-order WENO schemes. 

Next, a new one-parameter family of fifth-order WENO schemes is developed, 
and then is used as the starting point in the development of a family of “energy- 
stable” WENO schemes. Great care is exercised in the developing the WENO 
schemes to ensure that design-order accuracy is achieved in the vicinity of 
smooth extrema. 


3 Fifth- and Sixth-order WENO Schemes 


Any conventional high-order WENO finite- difference scheme for the scalar 
one- dimensional wave equation (1) can be written in the following semi discrete 
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form: 


dllj fj + | fj — | 

— ^ H 2- 2- = 0 

dt Ax 


(13) 


For the fifth-order WENO scheme that is presented in reference [2], the nu- 
merical flux fj + 1 is computed as a convex combination of three third-order 
fluxes defined on the following three-point stencils: Sll = { x j-2, x j-i,Xj}, 
Sl = {xj-i, Xj, Xj + i}, and Sr = {xj,Xj + \,Xj +2 }- (See figure 1.) Note that this 
set of stencils is not symmetric with respect to the (j + |) point; thus, the 
fifth-order WENO scheme is biased in the upwind direction. A central WENO 
scheme can be constructed from the conventional fifth-order WENO scheme by 
including an additional downwind candidate stencil Srr = {xj + i,Xj + 2 ,Xj +3 }, 
so that the collection of all four stencils is symmetric with respect to point 
(j + (;)• 1 The WENO flux, constructed in this manner, is given by 

fj+± = W j+l/2fj+l/2 + W j+l/2fj+l/2 + W f+l/2fj+l/2 + W f+l/2ff+l/2i (^) 


where , 2 , r = {LL, L , R , RR} are third-order fluxes defined on these four 
stencils: 


1 f LL (uj+ 1 / 2 ) ^ 
f L (u j+ 1 / 2 ) 
f R (u j+ 1 / 2 ) 


^ 2 —7 11 
-1 5 
2 

lo 


0 

2 

5 -1 
11 -7 2 


f(uj- 2) 

/K-i) 

/K) 

/(«i+i) 

/K+: 2) 

/K+3) 


(15) 


and ta LZ/ , ta L , and are weight functions that are assigned to the 
four stencils Sll, Sl, Sr, and Srr, respectively. The terms w LL , w L , w R , 
and w RR in equation (14) are nonlinear weight functions. These have both 
preferred values that are derived from an underlying linear scheme as well as 

1 The above approach to constructing central WENO schemes is proposed in refer- 
ence [6]. In general, central high-order WENO schemes that are built in this manner, 
are unstable when unresolved features or strong discontinuities exist in the compu- 
tational domain. The generalized energy-stable methodology presented in the next 
section guarantees the stability of the even order approximations while maintaining 
their nonoscillatory properties. 
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solution-dependent components. The preferred values are given by 


d LL = ■ d L = £ - 3p ; d R = ^ + 3p ; d RR = p , (16) 


where p is a parameter. The convergence rate of the scheme [eqs. (13)-(16)] 
with the preferred values w 6 J = d^ and r = {LL, L, R, RR} is greater than 
or equal to 5 for all values of the parameter p in equation (16). For the specific 
value p c = j, the fifth-order term vanishes and the convergence rate increases 
to 6 . The classical fifth-order upwind-biased WENO scheme of Jiang and Shu 
is obtained for p = 0 . 

The weight functions w LL , w L , w R , and w RR needed in equation (14) are 
given by 




(17) 


where 

a r 


d (r) 



t p \ 
e + /JM ) 


r = { LL , L , R , RR}, 


(18) 


and e < 0(Ax 2 ). The functions j3^ are the classical smoothness indicators. 
[See equation (59) for the general expression.] For fifth-order schemes, they 
are given by 

P LL = y§ (uj- 2 - 2 Uj_i + uj ) 2 + 1 (uj—2 ~ 4-Uj-i + 3 Ujf 

P L = if {Uj- 1 - 2u j + U 3+ if + \ (Uj- 1 - u j+1 f _ 

P R = y§ ( u j - 2u j+i + u j+zf + \ ( 3w i _ 4 m 1+i + u 3+ s) 2 

pRR = i| ( Uj+1 _ 2u i+2 + u i+3 ) 2 + 1 (— 5«j +1 + 8u i+2 - 3u i+3 ) 2 . 


The t p also varies with discretization order; the expression for 75 is given by 

t _ | (~fj- 2 + 5/j - 1 - 10/j + 10/j+i - 5/ j+2 + / j+3 ) 2 , for p ± 0 
I (fj-2 ~ 4/j-i + 6 fj - 4/ i+ i + / i+2 ) 2 , for P = 0 


Expressions for the fourth-, seventh-, and eighth-order WENO schemes appear 
in the appendix. 



The WENO mechanics expressed in equations (17)-(20) represents a signif- 
icant departure from the mechanics used in the original algorithm by Jiang 
and Shu [2], For example, equation (18) replaces the expression 


d (r) 

_ (6 + /JM) 2 


( 21 ) 


that is used in the weight functions given in reference [2], Furthermore, e < 
0( Ax 2 ) replaces e = O(10 -6 ) in [2], Equations (17)-(20) more closely resemble 
the weight and smoothness indicators proposed by Borges et al. in [ 8 ] There 
are differences, however, in how the parameters r p and e are chosen; differences 
motivated by the following observations. 

Borges et al.[ 8 ] proposes the following function f 5 : 

f 5 = \0 LL -j3 R \ ( 22 ) 


and a fixed value for the parameter e. Although the above weights and smooth- 
ness indicators [eqs. (17)— (19)] with the value of 75 given by equation (22) sig- 
nificantly outperform the conventional weights of Jiang and Shu [2] for both 
continuous and discontinuous solutions, three remaining shortcomings include: 

(1) The value of 75 given by equation (22) is not smooth because an absolute 
value function is used, which may reduce accuracy at points where (5 LL — 
/3 R changes sign. 

(2) This approach currently does not generalize to WENO schemes with 
design orders other than 5. For example, neither Tq = \/3 LL — j3 R | nor 
Tg = | j3 LL — (3 RR \ provides the design order of accuracy for the central 
sixth-order WENO scheme. 

(3) Near critical points where /', f", and f" vanish simultaneously, the mod- 
ified weights [eqs. (17)— (19) , and (22)] fail to provide the design order of 
convergence for the fifth-order WENO scheme if no constraints are im- 
posed on the parameter e other than e > 0. An example that demonstrates 
this property is presented in Section 6.1. 

Selecting r 5 in equation (20) to be a quadratic function of the fifth-degree 
undivided difference defined on the entire six-point stencil, circumvents the 
first, and second drawbacks encountered when using the scheme suggested in 
reference [8] . Indeed, 75 is a C°° function in its arguments, and can readily be 
generalized to WENO schemes of order p by choosing t p to be the /Ah- degree 
undivided difference that is defined on the entire (p + l)-point stencil. In 
contrast to the 75 found in reference [8], which is of order 0(Aa; 5 ) for smooth 
solutions, the proposed function r 5 as given by equation (20) is of order 0(Aa; 8 ) 
and 0(Aa: 10 ) for ip = 0 and tp ^ 0 , respectively; thus much faster convergence 
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of the fifth- and sixth-order WENO schemes to the corresponding underlying 
linear schemes is achieved. 2 

A remedy for the third shortcoming requires an additional modification to 
the approach proposed in reference [8]. Both the fifth- and sixth-order WENO 
schemes with the new weights that are given by equations (17-20) are design- 
order accurate for smooth solutions, including points at which the first- and 
second-order derivatives of the solution vanish simultaneously. However, if 
all derivatives up to fourth are equal to zero, then the fifth- and sixth-order 
WENO schemes locally become only third-order accurate. To fully resolve this 
issue, the constraint e < 0(Ax 2 ) is imposed herein. 

Equations (13)-(20) describe a new family of fifth-order WENO schemes. The 
primary difference between existing WENO schemes and this new family, is 
in the stencil biasing mechanics described by equations (17)-(20). Detailed 
theoretical justifications for the choice of r p and e used in equations (17)-(20) 
are presented in sections 5.2, 5.3. There, and again in the results section, it is 
shown that these parameters provide fast convergence of the new WENO and 
ESWENO schemes to the corresponding underlying linear schemes for smooth 
solutions, and deliver improved shock-capturing capabilities near unresolved 
features. 


4 A General Approach to Constructing High-Order Energy-Stable 
Schemes 


Algorithm (1) transforms an existing WENO scheme into an ESWENO scheme 
that is characterized by the following four properties: 1) a bounded energy 
estimate for arbitrary nonsmooth initial data, 2) conservation, 3) design order 
accuracy for sufficiently smooth data, and 4) discontinuity (shock) capturing 
capabilities that are similar to those of the base WENO scheme. 

By construction, algorithm (1) is applicable to 1-D periodic schemes of any 
order and produces a modified scheme that automatically satisfies the first and 
second properties: bounded energy estimate and conservation (See reference 
[1] for details.) Likewise, as shown in the next section, the additional terms 
added in the ESWENO formulation do not degrade the formal accuracy of 
the original WENO discretization (property three). However, the degree to 


2 The conventional fifth-order WENO scheme that corresponds to = 0, does not 
include the downwind stencil Srr] therefore, the fourth-degree undivided difference 
built on the available five-point stencil, is included in equation (20). This approach 
to handle the narrow stencil encountered with ip = 0, also generalizes to other orders 
of accuracy. 
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Algorithm 1 ( Transformation from WENO to ESWENO) 


1: Express the base WENO derivative operator as matrix D . 

2: Decompose D into symmetric and skew- symmetric components as D = D s k + D sym 

3: Add an artificial dissipation operator D a( i such that the modified symmetric matrix 

D S ym + D a d is positive semidefinite. 

■ T 

3a: Form the decomposition D sym = J2i=o[Di t ]Ai[Di' 1 } , 

2s + 1 is the bandwidth of the matrix D. 

The existence of this decomposition is established in the appendix. 

3b: Modify the diagonal terms of A* such that they are smoothly positive. 

One such approach is 



where Si, 0 < i < s are small positive constants that may depend on Ax. 


3c: Form the symmetric, positive semidefinite matrix 

D ad = P- 1 J2U 0 [D 1 i }A i [D 1 f (24) 

4: Form the energy-stable operator 

D = D+D a d (25) 


which the two formulations differ for unresolved data is not clear. Thus, the 
properties of the new ESWENO scheme must be tested to ensure it retains 
the desirable properties of the original formulation. 


ESWENO Schemes of Fifth- and Sixth-order 


We now apply algorithm (1) to the one-parameter family of fifth-order WENO 
schemes [eqs. (17)-(20)] for the scalar 1-D wave equation (1) with a > 0. 
Combining equation (15) with equation (14) and substituting the resulting 
WENO flux fj + 1 into equation (13) produces a stencil for the jth grid point 
of the form 


D 


5 

j.fc 


/ 0 



n ™f+ L l/2 

0 

0 

° \ 

0 

0 

-H+l/2 

5W j+l/2 

2w f+l/2 

0 

0 

0 

0 

0 

2w f+l/2 

5 “f+l/2 

~ lw f+i/2 

0 

0 

0 

0 

0 

llw RR 
j+ 1/2 

— 7 

(W j + 1/2 

2 

3 + 1/2 


7w i-i/a 

- llu fdi/2 

0 

0 

0 

0 

0 

1W f-l/2 


1/2 

0 

0 

0 

0 

0 

- 2 ™f-l/2 

- S -f-l/2 

lw ?-i/a 

0 

0 

V ° 

0 

0 

~ llw f R 1 / 2 

7w f- 1/2 

— 2 

3- 1/2 

0 / 


( f( u 3- 3) \ 
/(“j- 2) 

1 ) 

/(“j ) 

f( u J+ 2 ) 

\ /Oj+ 3) / 


(26) 
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with k on the interval j — 3 < k < j + 3. An explicit expression for the 
differentiation matrix D 5 follows immediately from equation (26). 


The derivative matrix D 5 is now decomposed into symmetric and skew-symmetric 
parts as 


D b 


D 


5 

skew 


+ D 


5 

sym * 


As with the third-order case (ref. [1]), the skew-symmetric component of D 5 
and the norm, take the forms 


n 5 = P -1 

skew 


Qz 


Qi 5 + Q§ — o 


P = A xl 


while the matrix D b is expressed as 


D 


sym 


= P 


1 ( D\ Al [£>?f + D\ Ag \D]} t + D\ \\ [£>;] t + £>» Ag [£>i>] T ) .(27) 


The matrices A) are diagonal with expressions for the jth element defined by 


(■^3 )j,j ~ 6 


,LL 


,RR 


W j+ 5/2 W j+l/2 


(28) 


( A 2 )j,j ~ 


12 


+«2 


4 «2 + <+3/2 

-<+ 1/2 + 4<? 1/2 - <+ 1/2 


(29) 


= £ 


+<-1/2 -<+1/2 


3<+l/2 “ 5 <i 3/2 + 2 «2 

+<+1/2 -<+3/2 

R 


+<-3/2 + 5 <-l/2 


_Qij « j RR 

6W j+ 1/2 


(30) 


wy = 


~ W j- 1/2 “ “f-1/2 “ “f-1 

,R 


3~ 1/2 

+«2 + <+ 1/2 + <+ 1/2 + <+ 1/2 . 


= 0 


(31) 


Because J2r w - r ' > — Aq is always equal to zero, and is not included in 
The sign of the diagonal terms (A* 5 ) - i = 1,2,3, could be either positive 
or negative; thus, the conventional fifth- and sixth-order WENO schemes may 
become locally unstable. If we define (A f) ■ z = 1,2, 3, to be smoothly positive 

= 2^ + ^ ) 

then the additional artihcial dissipation operator becomes D b d = P X)® =1 [-Di®] (A®) [Di] , 
and the resulting energy-stable scheme is obtained by adding the additional 
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dissipation term to the original WENO scheme. That is, 


D 5 = D s + Dlt . 


( 32 ) 


By construction, D 5 satisfies all of the conditions of Theorem 1, thereby pro- 
viding stability of the fifth- and sixth-order ESWENO schemes that are defined 
by equation (32). 


5 Consistency Analysis 

5. 1 Necessary and Sufficient Conditions for Consistency of WENO Schemes 

We now derive the necessary and sufficient conditions for the weight functions 
wf') for a family of pth-order WENO schemes to attain the design order of 
accuracy. Similar conditions have been obtained for the conventional fifth- 
order WENO scheme in reference [7]. With the same approach discussed in 
section 3 for p = 5, a one-parameter family of pth-order WENO fluxes can be 
constructed by using a convex combination of (s + 1) fluxes f (r ^ of stli order 
as 


fj± 1/2 = Z>j±i/ 2 /i±i/2» 

r 


(33) 


with 

fj±i /2 = K x 3 ± 1/2) + Hc^Aar' + 0{ Ax p+1 ), 

l=S 


(34) 


where h(x) is the numerical flux function that is implicitly defined as 


!{x) 



(35) 


and c 


(r) 

l 


are constants that do not depend on Ax. 


The corresponding pth-order WENO operator that approximates the first- 
order spatial derivative is given by 


[D p f]. 


f 7+1/2 ~ fj- 1/2 

Ax 



(r) Ar) 
j+l/2Jj+l/2 


— W 


(r) 

3 - 1/2 



Ax 


(36) 
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where [•],,■ is a jth component of a vector, the index r in equation (36) sweeps 
over all (s + 1) stencils, and is a nonlinear weight function that is assigned 
to the corresponding s-point stencil S r . For sufficiently smooth solutions, the 
weights wy> approach their preferred values cr r ), so that the WENO operator 
converges to the target linear operator _D Tai 's et as 


Target £• 

1 3 


r Target r Target 

Jj+l/2 Jj— 1/2 


Ax 


E (d^f. 


(r) 

3 + 1/2 


Ax 



(37) 


where d ^ is a one-parameter family of constants chosen to ensure pth-order 
convergence of the target operator to the exact value of the first-order deriva- 
tive at Xj. That is, 


Target £ 

= 

- 

3 dx 


+ 0{ Ax p ). 


(38) 


The coefficients d h) for one-parameter families of linear target schemes up to 
eighth order are given in the appendix. All target linear schemes in the family 
are (2s — l)th-order accurate, except one central scheme, which is (2s)th-order 
accurate. 


Subtracting equation (37) from equation (36) and using equation (34), we have 


[D p f 


J) Target f 

» 


E 




J j+ 1/2 ^ J — 1/2 


1 - 1/2 


J 3 


Ax 

( r ) 


S7 E [(Wj-+i /2 - d^)h (r \x j+ 1 / 2 ) - (w]E/ 2 - d (r) )/i (r) (x i _ 1 / 2 ) 


(39) 


+ £ E Aa/- 1 cJ r) [(fo+ 1/2 - S r) ) - (fo_’i /2 - + 0(Ax p ) 


J-lJO 


,(r) 


(r) 


From equation (39) it immediately follows that to retain pth-order accuracy, 
the weights of the WENO operator D p should satisfy the following necessary 
and sufficient conditions: 


E [w% 2 ~ d^\ = 0( Ax p+1 ) 


EcW [(«$ 1/2 - dV) - «_ j 1/2 - dW) J = 0 ( At ^- s+1 ) 


(r) 


(40) 



/ ( r ) 
K+l/2 


d^) — (w 


( r ) 

3 - 1/2 



O(Ax) 


Here, we use the following properties of h(x) and d*' : f\x) 
and E d^ = 1. 

r 


h(Xj+l/2)-HXj-l/2) 

Ax 
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By construction, the weights [eq. (17)] are normalized such that Y^r w - r ' > — 1; 
thus, the first constraint in equation (40) is satisfied identically. To simplify the 
analysis, especially for f{x) with an arbitrary number of vanishing derivatives, 
we hereafter use the following sufficient condition on wffi for the conventional 
WENO scheme to attain pth-order accuracy: 

w (r) _ d(r) = 0(Aa; p_s+1 ) . (41) 


This constraint is a direct consequence of the necessary and sufficient condi- 
tions [eq. (40)]. 


5.2 Sufficient Conditions for Consistency of ESWENO Schemes 


In this section, we show that the conditions [eq. (41)] and the constraints on 
Si in equation (23): 

5 t = Q( Ax p ~ 2i+1 ) , 1 < i < s , (42) 


guarantee that the energy-stable modifications of the conventional pth-order 
WENO scheme (see section 4) preserve the design order of the original scheme. 
As follows from equation (25), the ESWENO operator consists of two terms: 
one is the original WENO operator and the other is the additional artificial 
dissipation operator D ad that is given by equation (24). As shown in the 
previous section, the WENO operator is pth-order accurate if equation (41) 
holds. Hence, we need only show that D ad f = 0( Ax p ). 


Let us prove this conjecture for the one-parameter family of fifth-order ESWENO 
schemes that is presented in Section 4.1. Note that the term P^ 1 DfA 0 [Df] 7 
in equation (27) is identically equal to zero because of the normalization 
ffr — 1; it need not be included in D ad . Thus, the additional artificial 
dissipation term is given by 


D ad f = p- 1 J D 1 A 1 [ J D 1 ] r f + p-'D^Dlff + (43) 


Ai 



* = 1,2,3 


(44) 


where A,; is a jth diagonal element of the matrix A*. For simplicity, we have 
omitted the superscript 5 and the subscript (j,j) in this section. 

First, we evaluate the terms Ai, A 2 , and A 3 that are defined by equations (28- 
30). To simplify the derivation, the sufficient conditions [eq. (41)] for the one- 


15 



parameter family of fifth-order ESWENO schemes is rewritten in the following 
form: 


tt ,M _ S r > = (tp -- - i) 0(Ai 3 ) + 0( Ax 4 ) . 


(45) 


In equation (45), the stencil width s of each reconstruction polynomial is equal 
to 3, and the order p of this family of schemes is equal to 5 if (p ^ or 6 if 
p> = As follows from equation (45), the stiffer constraint on the weights 
should be imposed to obtain sixth-order accuracy. 

Replacing w ^ in A 2 with the corresponding preferred values d ^ and using 
equations (16) and (29), for any value of the parameter ip yields 


1 

12 


d LL - 4 d LL + d L -d R + 4 d RR - d RR 


= 0 . 


(46) 


Subtracting equation (46) from equation (29) and taking into account equation 
(45) yields 


Ao — T?7 


(wf+3/2 - d ' LL ) - 4 K>( 5/2 - d ' LL ) + K L +3/2 - dL 
-{wf+ 1/2 - d R ) + 4:(w rr 1/2 - d RR ) - (w RR lf 2 - d RR ) 
= (<p - 0{ Ax 3 ) + 0{ Ax 4 ) . 


(47) 


By comparing equations (29) and (30), one can see that Ai is at least one order 
higher than that of A 2 because of an additional cancellation that occurs within 
each group of terms that is associated with the same stencil. For example, 
expanding all of the terms w iL in equation (30) about Xj yields 


3 w 


LL 

I+1/2 


— 5 w 


LL 

3+ 3/2 


+ 2 w 


LL 

3+ 5/2 


= -2 


dw 


LL 


dx 


Ax + 0( Ax 2 ) 


(48) 


which gives an extra factor 0( Ax) compared with the corresponding terms 
w^ L in equation (30): — Aw RR 5 , 2 = 0(1). The same conclusion can be 

drawn for the other groups of terms that are associated with the L , R, and 
RR stencils in equation (30), which leads to 

A, = (v - T) 0(Ax 4 ) + 0(Ai 5 ) . (49) 
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Applying the same procedure to A 3 yields 


A3 


+ 


i ( d LL - d RR ) 
l(h5 - 2{ P) + (<d - 2h 


LL jLL\ f„..RR 

j+ 5/2 a ) \ W j+ 1/2 

) 0(Ai 3 ) + 0(Ai- 4 ) . 



( 50 ) 


To evaluate each term in equation (43), we consider three cases: 1) |A*| Si > 

0, 2) A i = 0(5i), and 3) |A,| <C <5*, i — 1, 2, 3. If we assume that |Aj| 3> > 0, 

then equation (44) can be expanded as follows: 


- _ | A*| + A i Sf 


(51) 


which yields 


A,- , if A,- > 0 


Ai 


Si , 

5 2 

l 


4|Ai| > 


if Ai = 0 , 
if Ai < 0 


(52) 


where the higher order terms have been omitted. Because equation (52) has 

been derived under the assumption that |A*| Si > 0, we can immediately 

s 2 

conclude that |A;| 3> Si S> Therefore, we only need consider that 

A; = A* , (53) 


which provides the lowest order of convergence for the term D\hi[D\ ] 7 f. If we 
substitute equation (53) in equation (43) and use equations (47, 49, 50), then 
the additional ESWENO dissipation term becomes 


Dadf — 


(v9- A)0(Aa; 3 ) + 0(Aa: 4 )] 

(<p - 0( Ax 2 ) + 0(Aa; 3 )] Dl\D\ 


(54) 


AO 


—2 <p 


6A# 


+ ~~ ^j) 0(Ax 2 ) + 0(Ax 3 ) D\[D\ 


■)3lTx 


If we take into account that D\[D\] T i = 0( Ax 21 ), then D ac ii can be recast as 
D ad f = (v - 0( Ax 5 ) + 0( Ax 6 ) . (55) 

From equation (55) we see, that the additional dissipation term is at least 
fifth-order accurate for all values of the parameter p . The fifth-order term 
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vanishes for ip = W for this value the order increase by one to sixth-order. 


The second case can be considered in a similar manner. Substituting A* = 
0(Si), i = 1,2,3, in equation (43) yields 


D ad f = + + 

= 0(8 1 Ax) + 0(S 2 Ax 3 ) + 0(8 3 Ax 5 ) . 


(56) 


To guarantee that the ESWENO dissipation term is pth-order accurate, the 
following constraints should be imposed on 8f 

<$i = 0(Ax p ~ 1 ) 

8 2 = 0(Ax p ~ 3 ) (57) 

8 3 = 0(Ax p ~ 5 ) , 

where p is equal to 5 for tp ^ ^ or 6 for ip = Note that the constraints [eq. 
(57)] are fully consistent with those that are given by equation (42). The third 
case, |Aj| -C Si, is similar to the second one and results in the same constraints 
[eq. (57)] on 8p therefore the third case is not presented here. 

Remark 5 Note that <5j, i = 1,2, 3, are user-defined parameters; therefore, 
the conditions [eq. (42)] can always be met. 

Remark 6 Although only fifth- and sixth-order ESWENO schemes have been 
analyzed in this section, the same procedure is directly applicable to the ad- 
ditional ESWENO schemes presented in the appendix. Thus, we can conclude 
that if equations (41) and (42) hold, then all the ESWENO schemes that are 
considered in this paper are design-order accurate. 


5.3 Consistency of the ESWENO scheme with New Weights 


The new weight functions for the one-parameter family of pth-order WENO 
and ESWENO schemes are given by 


( r ) 

■ur ; = 


a r 

E~a? 

1 


a r 


d {r) 



t p \ 
e + (30 ) ’ 


(58) 
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where (3 ^ are the classical smoothness indicators: 


/3« = gA^-‘ (59) 

x i~h 


q r (x) is an (s — l)th-degree reconstruction polynomial that is dehned on a 
stencil S r , and e is a small positive parameter that can depend on Ax. In 
equation (58), t p is dehned by 


(V Xj— S j r \, 1 . . 

■,Xj+ s >f , for tp ± 0 

(60) 

(V Xj— s -\-\ , . . 

■ ,Xj +s - 1 >) 2 , for tp = 0 

(61) 


where V < Xj- s+ i, . . . , Xj +S > is the pth-degree, undivided difference. Note 
that for the original WENO schemes of Jiang and Shu [2], which correspond 
to tp = 0, the entire stencil includes only (2s — 1) points; therefore, the highest 
degree of undivided differences that can be constructed on this stencil is 2s — 2, 
rather than 2s — 1, which is used for the other schemes in this one-parameter 
family. In particular, w^ r \ j3^ r \ and 75 for the fifth- and sixth-order WENO 
and ESWENO schemes are given by equations (17)-(20). Another scheme 
that requires special consideration is the central ESWENO scheme, which is 
obtained by setting tp = tp c , so that it is one order higher than that of the 
other schemes in the family. 

First, we present a truncation error analysis for the entire one-parameter fam- 
ily of pth-order ESWENO schemes, except for the two schemes that correspond 
to tp = 0 and tp = tp c . We demonstrate that the new weights that are defined 
by equations (58)-(60) satisfy the sufficient condition [eq. (41)] for smooth 
solutions with any number of vanishing derivatives if the following constraint 
is imposed on the parameter e in equation (58): 

e > 0{Ax p+s ~ 1 ) > 0, (62) 


where p is a design order of the scheme and (s — 1) is a degree of the corre- 
sponding reconstruction polynomials. 

If we assume that all of the required derivatives are continuous and use the 
properties of the Newton divided differences, then the Taylor series expansions 
of and t p ( tp ^ 0) at x 3 are given by 

/?( r ) = f' 2 Ax 2 + 0(Ax 3 ) , 


Tp = (f^ p) ) 2 Ax 2p + 0(Ax 2p+1 ) , 
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(63) 

(64) 



where f and /A> are the first- and pth-order derivatives of / at Xj. For ex- 
ample, for the family of fifth-order WENO schemes, these expansions are 


f3 LL = f 2 Ax 2 + ( A /" 2 - Iff") Ax 4 + (-f + iff"") Aar 5 + 0(Aa; 6 ) 

/ 3 L = f 2 Ax 2 + (A/" 2 + iff") Ax 4 + 0( Ax 6 ) 

P R = f 2 Ax 2 + (A /" 2 - iff") Ax 4 + ( A /"/'" - Iff"") Ax 5 + 0(Ax 6 ) 

(j RR = f 2 Ax 2 + ( A /" 2 - Aff") Aa;4 + (Af'f" _ 5f'f"") Ax 5 + 0( Ax 6 ) 

r 5 = (/ (5 ^) Aa; 10 + 0(Ax u ) , for p f 0 . (66) 


We first consider a case with f{xj ) f 0. Substituting equations (63) and (64) 
in equation (58) and accounting for equation (62) yields 

= 0(Ax^~ 2 ) (l - 0(A^ + >- 3 )) , (67) 

which leads to 

w (r) = d {r) + 0(Aa; 2p " 1 ) . (68) 


Note that the order of convergence of w to its preferred value is 2p — 1 rather 
than 2p — 2. The main reason for such “superconvergence” is the additional 
cancellation that occurs because the leading truncation error terms of all of 
the smoothness indicators are identical to each other if f f 0, as can be 
seen in equation (63). Equations (67) and (68) are valid only for p > 3 and 
are not applicable to the third-order WENO and ESWENO schemes that 
correspond to ip = 0. The detailed analysis of these third-order schemes (ip = 
0) is presented in reference [1]. By comparing equation (68) with equation 
(41), we can immediately conclude that the new weights satisfy the sufficient 
condition [eq. (41)], ensuring that both the WENO and ESWENO schemes are 
design-order accurate. Furthermore, for schemes of order 3 (ip f 0) or higher, 
the weights converge to their preferred values at a rate that is significantly 
faster than that given by the sufficient condition (41). The result is a much 
faster rate of convergence for both the ESWENO and WENO schemes to the 
corresponding target linear schemes, even on coarse and moderate grids. 

The next issue that is addressed is the convergence of the ESWENO schemes 
near the critical points at which the first-order and higher order derivatives of 
the flux approach zero. Let x c be a critical point at which the flux function is 
sufficiently smooth and at which its derivatives of order up to n V( )t\i are equal 
to zero. That is, 

f(x c ) = ... = f M (x c ) = 0 , f n ^ +1 \x c ) f 0 . 
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In contrast to the previous case for which f ^ 0, the leading truncation error 
terms of the smoothness indicators at the critical point are not equal to each 
other; thus, no additional cancellation occurs. For any number of vanishing 
derivatives, the following inequalities always hold: 

e > 0(Ax p+s - v ) > 0{ Ax 2p ) > t p , 


which yields 


e + /3W 


< 1 • 


By using equation (69), the weights can be recast as 

d (r) + d (r) ^ 


w {r) = 


l + E 

i 


= <jM + O 


Tr, 


:+A !) 


| jj( r ) I ’ 


(69) 


(70) 


where we use = 1- F° r an y number of vanishing derivatives, we have 


e + p<r) 



0(Ax 2 p) 

0(AxP+ s ~ 1 ) 


0{ Ax p - S+1 ) . 


(71) 


From equation (71) we see that u/ r ) converges to S’'' 1 at the rate of 0( Ax p s+1 ) 
or higher and satishes the sufficient condition (41). 

As mentioned at the beginning of this section, the two schemes that correspond 
to Lp = 0 and Lp = Lp c require special consideration. Using the same procedure 
that is outlined above, we can easily show that if the parameter e satishes the 
constraints 

e > 0(Aa; 3s-4 ) , for <p = 0 

(72) 

e > 0(Aa; 3s-3 ) , for p = p c , 


then the corresponding ESWENO schemes are design-order accurate regard- 
less of the number of vanishing derivatives of the solution. The constraints [eq. 
(72)] are derived using the following relations between the order of the scheme 
and the degree of the reconstruction polynomials: 


p = 2s — 1, 
P = 2s, 
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for p = 0 
for tp = p c . 



Also, note that if no constraint is imposed on e except that it must be strictly 
positive, then the order of convergence of the /All-order WENO and ESWENO 
schemes with the new weights may deteriorate from p to s. Indeed, for a 
sufficiently large number of vanishing derivatives n V( [ , (3^ and t p may become 
of the same order. For example, for the family of fifth-order schemes ^ 0), 
this type of degeneration occurs at n v d = 4. If f = f" = f" = f"" = 0, 
/ (5 ) 7^ 0, and e < 0(Aa; 10 ), which does not satisfy the condition given in 
equation (62), then equations (65) and (66) lead to 


r 5 0(Aa; 10 ) 

7TW = 0(Ax w ) + 0(Ax w ] 


= 0 ( 1 ) 


(73) 


As a result, the weights are of order 0(1) and the fifth-order WENO and 
ESWENO schemes locally degenerate to third order. 

Remark 7 In reference [8], the following modification of the weight functions 
given in equations (IT)— (19) , and (22) is proposed to recover the fifth-order 
rate of convergence if n v d < 2: 


i 


= d {r) 1 + 


A 

e + / 3 W 


(74) 


where m = 2. However, the modified weights [eq. (74)] still experience the 
same degeneration in accuracy near critical points for any choice of the pa- 
rameter m if n v d < 6. We demonstrate this with the following “thought” 
experiment. First, build a polynomial of degree 7 (or higher) such that at 
least 6 derivatives vanish at a single point. Next, connect this polynomial to a 
constant polynomial at the point where the derivatives of the first polymonial 
vanish. The solution to equation (1) is then a six times continuously differen- 
tiable function everywhere on the combined piecewise polynomials. Note that 
the underlying linear scheme is fifth-order accurate in this case. If we assume 
that the LL stencil is located completely in the constant part of the solution 
(i.e. fj _ 2 = fj _ i = fj ) while the other two stencils include points from the 
polynomial part of the solution, then we have /3 LL = 0, /3 L ^ 0, and f3 R ^ 0. 
As a result, 75 = /3 LL — /3 R = (3 R . If no constraint is imposed on e, then we 
can always choose e such that (3 R 3> e — > 0 on any given grid, which yields 

TO (R) _ d (R) = Q = O = 0(1) - O (jdj -> 0(1), 

This leads to w R = 0(1) and, consequently, to a loss of the design order 
of accuracy. Note that this problem persists regardless of the choice of m in 
equation (74). 
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Remark 8 The parameter e is user-defined, and therefore, the sufficient con- 
ditions [eqs. (62) and (72)] can always be satisfied. 

Remark 9 Equations (62) and (72) do not provide sharp estimates for the 
parameter e, which can be weakened if additional information regarding the 
number of vanishing derivatives is available a priori. Note, however, that the 
constraints [eqs. (62) and (72)] guarantee the design order of convergence of 
the corresponding WENO and ESWENO schemes for smooth solutions with 
an arbitrary number of vanishing derivatives. 


The last issue that we discuss in this section concerns discontinuous and unre- 
solved solutions. To successfully emulate the ENO strategy, a stencil for which 
the solution is discontinuous is eliminated from the approximation by effec- 
tively nullifying the corresponding weight that is associated with this stencil. 
Let a discontinuity be located inside a stencil S r , while the solution is smooth 
in all other stencils. We can easily verify that r p = 0(1), because t p involves 
all points of the entire stencil including those that contain the discontinuity. 
Therefore, the weight w *0 can be evaluated as 


fflW = 


d W 


1 + 


Q(b ' 

e+0(l)J 


0 ( 1 ) 


l^r 


0 ( 1 ) + 


0 ( 1 ) 


e+0(Ax 2 ) 


= 0(1) [e + 0(Aa; 2 ]) (75) 


Based on the above equation, to reduce the influence of e on the solution near 
the discontinuity, the parameter e should satisfy the following constraint: 

e < 0(Ax 2 ) . (76) 


Indeed, if equation (76) is met, then w ^ is of the order 0(Ax 2 ) and has the 
same order of magnitude as it would have if e = 0. Hence, the parameter e 
must be bounded not only from below by equations (62) and (72), but also 
from above by equation (76). 

Another consideration that can help us optimally select the parameter e is the 
manner in which the WENO and ESWENO schemes handle small-amplitude 
oscillations. Consider a solution that contains small-amplitude spurious oscil- 
lations 


f = f, + Sf, 


(77) 


where f s is a smooth component and Sf is a nonsmooth, high-frequency com- 
ponent of the solution. By substituting equation (77) in equation (59) and 
assuming that the contribution of the smooth component is negligibly small, 


23 



we have 


/3 (r) = 0(6 f 2 ) . 

As follows from equation (58), if e 3> 0(6 f 2 ), then the denominator e + j3 ^ 
is dominated by e, and these spurious oscillations cannot be detected by the 
ESWENO dissipation mechanism. However, if e < 0(6f 2 ), then the weights 
begin to deviate from their preferred values, which increases dissipation in 
those stencils containing spurious oscillations. 


From these considerations it follows that the lower bound of the constraints on 
e given by equations (62), (72), and (76) should be used (1) to obtain the design 
order of convergence at the smooth parts of the solutions, (2) to effectively 
damp high-frequency spurious oscillations, and (3) to provide good shock- 
capturing capabilities near strong discontinuities. Therefore, in onr numerical 
experiments, we use 


e 


0(Aa; 3s-4 ) , for (p = 0 
< 0( Ax 3s ~ 3 ) , for (p = (p c 
0(Ax 3s ~ 2 ) , otherwise 


(78) 


where s — 1 is the degree of the reconstruction polynomials. With the above 
selection of e, all spurious oscillations of amplitude 0(e 1 ^ 2 ) and higher are 
suppressed by the ESWENO dissipation mechanism. The same conclusions 
can be drawn for the WENO schemes with the new weights given by equations 
(58)— (61). 


6 Numerical Results 


We now assess the performance of the fourth-, fifth-, and sixth-order ESWENO 
schemes with the new weights and compare them with the conventional WENO 
counterparts. For all of the numerical experiments that are presented, the pa- 
rameters e of the ESWENO weight functions and the parameters <5* of the 
additional artificial dissipation operator are chosen based on equations (78) 
and (42) with two minor modifications. First, to preserve the scale invariance 
of the original WENO schemes, the physical grid spacing Ax in equations (42) 
and (78) is replaced with A£, where A£ = 1 / J and J is the total number of grid 
cells. Second, as follows from equation (58), the parameter e should be scaled 
consistently with /Wh In regions where the solution is smooth, for the third- 
and fourth-order WENO and ESWENO schemes, (3^ approximates f' 2 A£ 2 . 
For the fifth- and sixth-order schemes, /T r ) approximates the linear combi- 
nation of the same first-order derivative term and a second-order derivative 
term that is proportional to f" 2 . The same pattern persists for higher order 
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exact 



Fig. 2. Solutions obtained with fourth-order WENO and ESWENO schemes on 
101-point grid for linear wave equation with initial condition [eq. (81)] at t = 0.04. 

schemes as well. In the vicinity of discontinuities, f^ r ' 1 is of the order 0(/ 2 ). If 
we take into account the above considerations, for all test problems considered 
the parameter e is set to 


e 


C 


CA£ 3s_4 , for (p = 0 
< CA£ 3s ~ 3 , for (p = (p c 
C A£ 3s_2 , otherwise 



(79) 


where /o = /[wo(0] i s the initial flux, w 0 (£) is the initial condition, f{f is 
the (s — l)th-order derivative of f 0 with respect to £, s — 1 is the degree of 
the reconstruction polynomials, Ai is a set of points at which the solution is 
discontinuous, and || • || is a norm in which the solution is sought. The scaling 
factor C in equation (79) can easily be evaluated because it depends only on 
the initial condition for which the locations of all discontinuities are known 
a priori. Note that the parameter e is calculated once and the same value is 
used over the entire time interval of integration; thus the computational cost 
does not increase. 

In accordance with equation (23), the parameter S t should be scaled consis- 
tently with A,;. Because A* is a linear combination of the weights with each 
being of the order 0(1), the scaling factor of S t is set equal to one, which leads 
to 


Si = A£ p_2i+1 . 


( 80 ) 
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Equations (79) and (80) eliminate the ambiguity in determining the parame- 
ters e and Si for the ESWENO schemes; thus, the equations for e and Si are free 
of tuning parameters. Furthermore, equations (79) and (80) are fully consis- 
tent with the sufficient conditions [eqs. (42) and (78)] and allow the ESWENO 
schemes to be invariant when the spatial and time variables are scaled by the 
same factor. Note that the parameter e for the conventional WENO schemes 
of Jiang and Shu is set to 10 -6 as recommended in [2], In accordance with 
the recommendations of Borges et al. [8], the parameter e for the fifth-order 
WENO scheme with the weights given by equations (17) — (19) , and (22), which 
is referred as WENO-Z, is set to 10 -40 . 

The time derivative for all steady test problems is approximated by using a 
third-order total variation diminishing (TVD) Runge-Kutta method that is 
developed in reference [11], while unsteady problems are integrated by using a 
fourth-order low-storage Runge-Kutta method (ref. [12]). To reduce the fourth- 
order temporal error component and make it consistent with the spatial error 
of fifth- and sixth-order schemes, the time step in global grid refinement studies 
is reduced by a factor of 2 6//4 for each doubling of the number of grid points 
in space. The Courant-Friedrich-Levy (CFL) number has been set to 0.3 and 
0.6 for the steady and unsteady test problems, respectively. 


6. 1 Scalar Linear Wave Equation 


We begin by verifying that the new class of ESWENO schemes provide the 
design order of convergence for smooth problems, including local extrema. To 
check this property, we consider equation (1) with a = 1 and the following 
initial condition: 

Uq(x) = e - 30 °O-^) 2 , (81) 


where x c is 0.5. The computational domain for this test problem is set to 
0 < x < 1. Numerical solutions are calculated on a sequence of globally 
refined uniform grids and advanced in time up to t = 1, which corresponds to 
one period in time. 

First, we show that the conventional fourth-order WENO scheme is unsta- 
ble for this smooth problem; however, the corresponding ESWENO is stable, 
and its solution is in excellent agreement with the exact solution, as shown 
in Figure 2. To make the conventional fourth-order WENO scheme stable, 
a smoothness indicator that corresponds to the downwind stencil has been 
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Fig. 3. Lqo error norms obtained with the fourth-order linear, WENO, and 
ESWENO schemes for linear wave equation with initial condition given in equa- 
tion (81). 

modified as follows: 


0 R = 


{0 L ) k + (0 c ) k + {0 R ) k 


1 1 /* 


(82) 


where A; is a constant that is greater than 1. Note that /3 R is a C°° function 
of its arguments and approaches rna x(/3 L , /3 C , (3 R ) as k — > oo. In contrast to 
equation (82), a modified smoothness indicator /3 R = max(/3 L , j3 c , (3 R ) that is 
proposed in reference [6] is a nonsmooth function of /3 L , (5° and /3 R , which may 
lead to the degeneration of the design order of accuracy and is more prone 
to spurious oscillations near unresolved features and strong discontinuities. 
For k — > oo, the downwind smoothness indicator that is given by equation 
(82) prevents the corresponding weight w R from being larger than the lesser 
of the other two weight functions; thus, the stencil is biased in the upwind 
direction near unresolved features. In smooth regions, all three smoothness 
indicators are of the same order, and the weights approach their preferred 
values of w L = w c = |, and w R = which provides the design order of 
accuracy. For all of the numerical experiments that are presented herein, the 
parameter k in equation (82) is equal to 4. 

Figure 3 shows the L ^ error norms that are obtained with the fourth-order 
WENO and ESWENO schemes, and the corresponding underlying linear schemes. 
As shown in Figure 3, the fourth-order ESWENO scheme is significantly more 
accurate than the conventional fourth-order WENO scheme. The L <*, error 
that is obtained with the fourth-order ESWENO scheme is slightly higher 
than that obtained with the corresponding linear scheme on coarse meshes 
and reaches its theoretical limit starting at J — 201 grid points. The max- 
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Fig. 4. Loo error norms obtained with the fifth- (left) and sixth-order (right) linear, 
WENO, WENO-Z, and ESWENO schemes for linear wave equation with initial 
condition given in equation (82). 

imum error occurs at the peak of the Gaussian pulse indicating that the 
ESWENO scheme is design-order accurate at the smooth extrema. In con- 
trast to the ESWENO scheme, the conventional WENO scheme demonstrates 
only a third-order convergence rate, even on the finest mesh with J = 1601, 
and is two to three orders of magnitude less accurate on moderate and fine 
meshes. 

If the new weights [eqs. (58) and (59)] are used instead of their conventional 
counterparts, then the design order of convergence of the fourth-order central 
WENO scheme is recovered, as shown in figure 3. Moreover, the fourth-order 
WENO scheme with the new weights provides slightly better accuracy on 
coarse meshes than the corresponding ESWENO scheme. This result is not 
surprising because the ESWENO scheme has the additional dissipation term 
[eq. (24)], which guarantees the stability. In general, if a WENO scheme with 
the new weights is stable, then it is slightly less dissipative than the corre- 
sponding ESWENO counterpart. For all of the problems that are considered, 
however, the results obtained with the conventional WENO schemes with the 
weights given by equations (58) and (59) are practically indistinguishable from 
those of the corresponding ESWENO schemes; therefore these results are not 
presented hereafter. 

The Loo error norms that are obtained with the fifth-order WENO-Z scheme 
and the fifth- and sixth-order WENO and ESWENO schemes for the same 
test problem are depicted in figure 4. Similar to the fourth-order case, the 
conventional sixth-order central WENO scheme is unstable. This problem is 
avoided by using the same upwinding technique that is outlined earlier. For 
the sixth-order central WENO scheme, the modified smoothness indicator that 



corresponds to the most downwind stencil is given by 


fj RR = 


(/ 3 LL ) k + (/ 3 L ) k + (/ 3 R ) k + (/ 3 RR ) k 


l/k 


( 83 ) 


Both the hfth-order WENO-Z and fifth- and sixth-order ESWENO schemes 
are equal in accuracy to the corresponding underlying linear schemes, as shown 
in Figure 4. Although the fifth-order WENO scheme exhibits the design-order 
convergence rate, its L <*, error norm is nearly an order of magnitude larger 
than those of the corresponding fifth-order WENO-Z and ESWENO schemes. 
In contrast to that of the sixth-order ESWENO scheme, the L ^ error that is 
obtained with the conventional sixth-order WENO scheme shows only fifth- 
order convergence and is quite far from the theoretical limit (represented by 
the corresponding sixth-order central linear scheme). Note that the WENO-Z 
schemes of fourth- and sixth-order are currently unavailable in the published 
literature and therefore are not presented herein. 


5th-order WENO 
5th-order ESWENO 



Fig. 5. Time histories of rightmost eigenvalue of symmetric part of fifth-order WENO 
and ESWENO operators for Gaussian pulse problem. 

As we have proven in sections 2 and 4.1, for the linear convection equation (1) 
with piecewise continuous initial conditions, the family of ESWENO schemes 
is stable in the energy norm. The sufficient condition for stability is the neg- 
ative semidefiniteness of — \{D + D T ), where D is the ESWENO discrete op- 
erator. This property implies that all eigenvalues of the symmetric portion of 
the ESWENO operator are nonpositive. In contradistinction to the ESWENO 
scheme, the symmetric portion of the WENO operator may have positive 
eigenvalues. (See section 4.1.) Note that the same conclusion can be drawn for 
the fifth-order WENO-Z scheme (ref. [8]), whose derivative operator is identi- 
cal to that of the conventional WENO scheme with the modified weights [eqs. 
(IT) — (19) , and (22)] which vary in the interval from 0 to 1 near the unresolved 
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features. These properties of the WENO, WENO-Z, and ESWENO schemes 
are shown in figure 5, which gives the time histories of the rightmost eigen- 
value of the symmetric portion of the operators computed on a 201-point grid. 
The rightmost eigenvalue of the ESWENO operator — |(D + D T ) is equal to 
zero up to the order of the round-off error, while the symmetric part of the 
conventional fifth-order WENO and WENO-Z operators have positive eigen- 
values of order O(10 _1 ) and 0(10), respectively. (See figure 5.) This results 
from the fact that the symmetric portion of these WENO-type operators is 
not negative semidefinite, if unresolved features exist in the computational do- 
main. Note that the presence of the positive eigenvalues does not imply that 
the fifth-order WENO and WENO-Z schemes are globally unstable because 
these eigenvalues correspond to different grid points at different moments of 
time. 



Fig. 6. Lao error norms obtained with fifth-order linear, WENO, WENO-Z, and 
ESWENO schemes for linear wave equation with initial conditions given in equation 
(84) at t = 1. 


The superiority of the ESWENO schemes with the new weights compared 
with the conventional and modified WENO counterparts, becomes evident 
when a more challenging problem is considered. In the previous test problem, 
the initial condition [eq. (81)] has the critical point at which f — 0 and f" ^ 0 
(i.e., the number of vanishing derivatives is n v d = 1). As shown in section 5.3, 
the fifth-order WENO-Z scheme that was developed by Borges et al. [8] fails 
to recover the design order of convergence if n V( -j > 6 for an arbitrary choice 
of the parameter m in equation (74). We now demonstrate this property for 
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n v d > 3. Consider equation (1) with the following initial condition: 


u Q (x) = { 


z i& _ 14z i6 + 69 ,i4 _ 175z i2 + 2 59 -io 

— 231z 8 + 119^ 6 - 29^ 4 + 1, for |z| < 1 

0 otherwise 


( 84 ) 


where z = 5 (a; — 0.5) and 0 < x < 1. The above function is six times continu- 
ously differentiable and has three critical points: one at x — 0.5 with n vc i = 3 
(i.e., /'(0.5) = /"( 0.5) = /"'(0.5) = 0 and /""(0.5) ^ 0) and the remaining 
two at x — 0.3 and 0.7 with n V( i = 6. If we compare the L ^ error norms com- 
puted with the fifth-order linear, WENO, WENO-Z, and ESWENO schemes 
at t — 1, (see fig. 6), we see that the situation changes dramatically as com- 
pared with the previous test problem. As expected, the WENO-Z scheme fails 
to deliver fifth order convergence, while the ESWENO scheme is fifth-order 
accurate and provides practically the same error convergence as the underly- 
ing linear scheme. This result is not surprising because the parameter e for 
the WENO-Z scheme is set to 10 -40 , as suggested in reference [8]. As a result, 
at a point x = 0.5 + Ax the weights become of order 0(1). Thus, if we use 
equations [(65), (22), and (74)] with e = 0 and take into account that for the 
polynomial in equation (84), f = 0(Aa; 3 ), f" = 0( Ax 2 ), f" = 0( Ax), and 
f"" = 0(1) at x = 0.5 + Ax, then we have 


w (r) _ d (r) = q 


- |^ rr _ /TW | Aa; 5 - 

J' 2 Ax 2 + (ff' 2 - f’f'")Ax\ 


( 0(Aa: 8 )\ 

yO(Aa: 8 ) J 


0(1) .(85) 


As follows from the above estimate, the WENO-Z scheme fails to recover fifth 
order convergence regardless of the choice of the parameter m in equation 
(74). Our numerical results corroborate this conclusion. Figure 6 shows that 
increasing the exponent m in equation (74) results in an even larger dete- 
rioration in accuracy. In contrast to the WENO-Z scheme, the conventional 
fifth-order scheme of Jiang and Shu[2] converges at the design order and be- 
gins to approach the theoretical limit on the finest meshes. The main reason 
for such a behavior is the choice of the parameter e, which is set to 10” 6 . As 
the grid is refined, /r r - ) — > 0, the denominator in equation (21) is dominated 
by e, and the weights approach their preferred values. 


6.2 The 1-D Euler Equations 


In reference [1], the third-order ESWENO scheme is proved to be energy-stable 
for a system of linear hyperbolic equations with periodic boundary conditions. 
This result can be directly extended to the class of higher order ESWENO 
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schemes that are presented in this paper. For nonlinear conservation laws with 
nonperiodic boundary conditions, a similar proof for the high-order ESWENO 
schemes is not currently available. Nevertheless, we would like to test how the 
new schemes perform for the system of the quasi-l-D Euler equations, which 
are given by 


9U | <9F 
dt ' dx 


p 


pu 


pu 

pu 

, F = 

pu 2 + P 

0 

1 

pu 2 

E 


(. E + P)u 


(E + P)u 


( 86 ) 


P = (7-1) (£ + *£) , 


where A = A(x) is the cross-sectional area of a quasi-l-D nozzle and 7 = 1.4. 
As has been shown in reference [1], the ESWENO reconstruction must be 
implemented in local characteristic fields to guarantee the stability of the 
ESWENO scheme for systems of hyperbolic conservation laws. In the present 
analysis, both the WENO and ESWENO reconstructions are based on the 
Lax-Friedrichs flux splitting. See reference [1] for further details. 


3rd-order linear 5th-order linear 




Fig. 7. Lqo error norms obtained with the third- and fourth-order (left) and fifth- 
and sixth-order linear, WENO, and ESWENO schemes for the subsonic quasi-l-D 
nozzle problem. 

For the 1-D Euler equations, the preferred biasing of the stencil in the upwind 
direction, given by equations (82) and (83), is used for the central fourth- and 
sixth-order WENO and ESWENO schemes to suppress spurious oscillations 
near strong discontinuities. In the ESWENO case, these oscillations can also 
be eliminated by increasing the coefficient in front of the D\AiDj term in the 
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artificial dissipation operator. However, this approach is more dissipative and 
also reduces the maximum CFL number for which the scheme remains stable. 
For the test problems that are presented in this section, the results obtained 
with the fifth-order WENO-Z scheme are practically indistinguishable from 
those of the fifth-order ESWENO scheme; therefore the WENO-Z results are 
not presented hereafter. 

To verify that the new ESWENO schemes are design-order accurate for hyper- 
bolic systems, we consider the steady-state isentropic flow through a quasi- 
1-D nozzle with the following cross-sectional area: A(x) = 1 — 0.8a;(l — x), 

0 < x < 1, as a test problem. The inflow Mach number is set to 0.5 and 
the pressure at x — 1 is assumed to be equal to that at x — 0. Under these 
conditions, the flow is fully subsonic, and the solution is smooth. Global grid- 
refinement studies for the third-, fourth-, fifth-, and sixth-order WENO and 
ESWENO schemes are presented in figure 7. 

The Lqo error norms that are obtained with the third- and fourth-order ESWENO 
schemes exhibit the design order of convergence. On the finest mesh with J = 
801, the fourth-order ESWENO solution is approximately three and four or- 
ders of magnitude more accurate than those of the third-order ESWENO and 
WENO schemes, respectively. As shown in figure 7, the fifth-order ESWENO 
scheme provides the same accuracy as its underlying linear scheme, thereby 
exhibiting the perfect error convergence for this test problem. Although the 
conventional fifth-order WENO scheme exhibits the design order of conver- 
gence, its Lqo error norm is larger by a factor of 3 than that obtained with 
ESWENO counterpart. Similar to the fourth-order case, the central sixth- 
order ESWENO scheme converges at the design-order rate on both coarse 
and fine meshes, thus providing significantly more accurate solutions than 
both fifth-order schemes. 

Remark 10 For this problem, the conventional central fourth- and sixth- 
order WENO schemes do not converge to a steady-state solution even with 
the upwinding mechanisms [eqs. (82) and (83)] turned on. The primary reason 
for this tendency is the presence of positive eigenvalues in the spectrum of the 
symmetric part of the conventional WENO operator. 

The next test problem is the steady transonic flow through a quasi- 1-D nozzle 
with the following cross-sectional area: 

A{x) = 1.398 + 0.347 tanh(0.8x — 4), 0 < x < 10. 

The Mach number at x — 0 is 1.5, and the outflow conditions have been 
chosen so that the shock is located at x — 5. Density profiles are calculated 
using the fifth-order WENO and fifth- and sixth-order ESWENO schemes on 
a 51-point grid and are compared with the exact solution in figure 8. For all 
schemes, the captured shock is smeared over two grid cells, and the numerical 
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Fig. 8. Comparison of fifth-order WENO and ESWENO schemes for steady tran- 
sonic flow through quasi-l-D nozzle. 

solutions are essentially nonoscillatory and agree quite well with the exact 
solution. Note that the fifth- and sixth-order ESWENO schemes converge to 
the steady-state solution an order of magnitude faster than the conventional 
fifth-order WENO scheme; this result again indicates the presence of unstable 
modes that are generated by the positive eigenvalues of the WENO dissipation 
operator. 




Fig. 9. Density profiles computed with the fourth-order (left) and sixth-order (right) 
WENO and ESWENO schemes on a uniform grid with 201 points for the Sod 
problem at t = 0.16. 


The last two problems considered are standard unsteady problems for testing 
shock-capturing schemes. The first problem is the Riemann problem with the 
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initial conditions proposed by Sod[13]: 


(P,U,P) 


(1,0,1) if -0.5 < x < 0 

(0.125,0,0.1) if 0< x <0.5 


The numerical solutions that are computed with both fourth- and sixth-order 
WENO and ESWENO schemes are essentially nonoscillatory, as one can see 
in figure 9. Note, however, that the fourth- and sixth-order ESWENO schemes 
provide better resolution near the contact discontinuity and the shock as com- 
pared with the conventional WENO schemes. 

We conclude this section by presenting the numerical results that are obtained 
with WENO and ESWENO schemes for the shock entropy- wave interaction 
problem. The solution of this benchmark problem contains both strong dis- 
continuities and smooth structures, and is well suited for testing high-order 
shock-capturing schemes. The governing equations are the time-dependent 
1-D Euler equations [eq. (86)] with G = 0, subject to the following initial 
conditions: 


{p,u,p) 


(3.857134,2.629369,10.33333), if - 5 < x < -4 

(1 + 0.2 sin 5x, 0, 1) , if - 4 < x < 5 . 


(87) 


The governing equations are integrated in time up to t — 1.8. The exact 
solution to this problem is not available. Therefore, a numerical solution that 
is obtained with the conventional fifth-order WENO scheme on a uniform grid 
with J = 4001 grid points is used as a reference solution. 

Numerical solutions that are computed with the fourth-, fifth-, and sixth- 
order WENO and ESWENO schemes on a 301-point grid at t — 1.8 are com- 
pared with the “exact” reference solution in figure 10. Both the WENO and 
ESWENO solutions are free of spurious oscillations. The fourth-order WENO 
scheme is the most dissipative among the considered schemes and has the 
worst resolution in a region just upstream of the moving shock. The fourth- 
order ESWENO scheme provides better resolution of the high-frequency oscil- 
lations behind the shock. However, the wave amplitudes are much lower than 
those that are computed with the fifth- and sixth-order schemes. Again, the 
fifth- and sixth-order ESWENO solutions are more accurate than those of the 
corresponding WENO schemes. Both ESWENO schemes resolve the smooth 
extrema quite well, and the solutions are essentially nonoscillatory near the 
captured shock. 
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exact exact 



Fig. 10. Density profiles computed with the fourth-order (top left), fifth-order (top 
right), and sixth-order (bottom) WENO and ESWENO schemes on uniform grid 
with 301 grid points for density sine wave/shock-interaction problem. 

7 Conclusions 


A systematic methodology is developed to construct Energy Stable weighted 
essentially nonoscillatory (ESWENO) finite- difference schemes of arbitrary or- 
der, and is used to construct ESWENO schemes based on the WENO schemes 
that are presented in references [2] and [6]. The new ESWENO schemes differ 
from the conventional WENO schemes by the addition of a nonlinear artificial 
dissipation term of special form. The additional term is design-order accu- 
rate for smooth solutions, including smooth extrema, and guarantees that the 
ESWENO scheme is stable in the /^-energy norm for both continuous and 
discontinuous solutions of hyperbolic systems; for the conventional WENO 
scheme, an energy estimate is not available. The distinctive feature of the new 
class of ESWENO schemes is that the spectrum of the symmetric part of the 
ESWENO operator is always located in the left half-plane, while the symmet- 
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ric part of the WENO operator has positive eigenvalues; thus, the conventional 
WENO schemes may become locally unstable. 

New weight functions are also developed, guided by newly derived sufficient 
conditions guaranteeing design order accuracy for arbitrary WENO schemes; 
sufficient conditions that are valid regardless of the number of vanishing so- 
lution derivatives. A truncation error analysis is used to optimize the weight 
functions, so that the new ESWENO schemes satisfy the sufficient condi- 
tions for accuracy and provide essentially nonoscillatory solutions for problems 
with strong discontinuities. Numerical experiments confirm that the high-order 
ESWENO schemes with the new weights are significantly more accurate than 
conventional WENO schemes of the same order, and demonstrate excellent 
shock-capturing capabilities. 


8 Appendix 

8. 1 Existence of Symmetric Decomposition 

The following lemma establishes the existence of the decomposition D sym = 
X^ = o[£V]A i{Df] . The lemma provides an algorithm to decompose a matrix 
of arbitrary bandwidth 2s + 1 into the sum of two matrices: one that can be 
expressed as one term of the desired symmetric form, and one of bandwidth 
2(s — 1) + 1. 

Lemma 11 Define E s to be an A- dimensional, symmetric, bidiagonal matrix 
with the subdiagonal elements e s+ i t i, e s+2 ,2, •••, ejv-i,jv-a-i, e N,N- s - The 
parameter s is the offset from the main diagonal. Further, define B s _i to be 
an A- dimensional, symmetric, banded matrix of bandwidth 2 (s — 1) + 1 with 
elements that are functions of the matrix E s . The matrices E s and B s _\ satisfy 
the following relationship: 

E s = + B,-! . 


with A s = diag[e s+ i } i,e s+2 , 2 r ■ ■ ,eN-i,N-s-i,e N ,N-s,O s T ] and 0 S is a zero 
vector of dimension s. 

Proof: A proof by induction is not included herein. One step of the algorithm 
can be verified immediately by inspection. 

A simple example illustrates Lemma 11. Define the 5x5 matrices E- 2 , D 15, 
and B\ such that E 2 = (— l) 2 [Di5 2 ]A 2 [Di 5 2 ] T + B\ and 
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e 2 


1 0 0 a 0 0 N 
0 0 0 6 0 
a 0 0 0 c 
0 6 0 0 0 
^0 0 C 0 0y 


D 


15 


^ 1 0 0 0 - 1 ^ 

-110 0 0 
0-1100 
0 0-110 
v 0 0 0 -1 lj 


a 2 


a 0 0 0 0^ 


( 

a 

—2a 

0 

0 

o N 

0 6 0 0 0 



—2a 

4a + 6 

— 2(a + 6) 

0 

0 

0 0 c 0 0 

; B x = 


0 

— 2(a + 6) 

a + 46 + c 

-2(6 + c) 

0 

0 0 0 0 0 



0 

0 

-2(6 + c) 

6 + 4c 

-2c 

0 0 0 0 0y 


l 

0 

0 

0 

-2c 

c ) 


Remark 12 Lemma 11 can be used recursively, beginning with the outermost 
s diagonal and working inward, to establish that D sym = Y^i=o[Di l ]Ai[Di] . 
Each step generates a new symmetric matrix with a bandwidth that is equal 
to 2 less than the previous value until only a diagonal matrix remains. Thus, 
the algorithm terminates after s + 1 steps with the desired decomposition. 

Remark 13 This recursive algorithm works for nonperiodic matrices in mul- 
tiple spatial dimensions. 


8.2 ESWENO Schemes of Third- and Fourth-Order 


Based on equation (13), /- + 1 defined for a third- or fourth-order WENO 
scheme is 

fj+k = W j+l/2fj+l/2 + W j ; +l/2fj+l/2 + W f+l/2fj+l/2 ( 88 ) 


where w L , w c , and w R are weight functions that are assigned to each respec- 
tive stencil. The second-order linear fluxes: /j+i/ 2 > defined in equation (88) for 
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r = {L, C , i?}, are 


f f L (uj+ 1/2) ^ 
f C (u j+ 1/2) 

V / R («j+l/2) ) 


<13 0 0 N 
0 110 
v 0 0 3 -1, 


1 f(Uj - < 
f(Uj) 
f(Uj+ 1 ) 

v/K+2); 


(89) 


The terms , ic*- 7 , and are the nonlinear weight functions that are given 
by equations (58) and (59) with 

T _ \ (“/i-i + 3 /j ~ 3 /j+i + fj+ 2) 2 ’ for V ^ 0 
[ (./) 1 - 2/j + /j+ 1 ) 2 , for v? = 0 

(See ref. [1] for further details). These weight functions have preferred values 
that are derived from underlying linear schemes, as well as solution-dependent 
components. The preferred values are given by the formulae 

d, L = | — p ; dp = | ; d R = p , (90) 


where p is a parameter. The convergence rate of equation (13), with the pre- 
ferred weight values that are given in equation (88), is equal to 3 for all values 
of the parameter p except p c — |, for which the convergence rate is 4. 


Combining equation (89) with equation (88) and substituting the resulting 
WENO flux fj + 1 into equation (13) produces a stencil for the jth grid point 
of the form 




' 0 

“<+1/2 

3 <+ i /2 

0 

0 N 


0 

0 

<+1/2 

<+1/2 

0 

1 

0 

0 

0 

3<+l/2 

-<+1/2 

2 Ax 

W J - 1/2 

~ 3 w )- M 2 

0 

0 

0 


0 

-<-1/2 

-<-1/2 

0 

0 


l 0 

0 

-3<- i /2 

<1/2 

0 ) 


OK- 2 )^ 

f(uj- 1 ) 
f(Uj ) 
/K+i) 
\f(Uj+ 2 ) ) 


(91) 


with k on the interval j — 2 le k < j + 2. 
Simplifying D 3 using the expression 

w c = 1 - w L - w R 
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yields expressions for a single row of the D 3 matrix of the following form 

T 


( 


n 3 


2Ax 


\ 


0 


W J- 1/2 


- 2w j- 1/2 - <+1/2 + <-1/2 - 1 


<-1/2 + 2 <+l/2 


2<-i/ 2 


w 


R 


3+ 1/2 


A 

j+1/2 


-wr, 1 /O + <Li/2 + 2^^ /9 + 1 


,R 

J j+ 1/2 


-<+1/2 


V 


) 


( 92 ) 


The derivative matrix Zb 3 is now decomposed into symmetric and skew-symmetric 
parts 

n 3 _ +)3 I n3 

skew 1 sym 

As with the fifth-order case, the skew-symmetric component and the norm of 
D 3 take the forms 

Dskew — p <3 ; Q3 + Ql = 0 ; p = ■ 

while the matrix D'l ym is expressed as 

D 3 „ m =P- l (D 1 2 A 3 2 \D 1 2 f +D 1 1 Afp 1 1 ] T + D 1 “Aj[D 1 °] T ) (93) 


where A 3 is a diagonal matrix with the expressions for the jth element defined 
by 




3,3 4 L 

1 r 


<+3/2 - <+1/2 


( A l)/4 = 4 r W 3+ 3/2 + W J+l/2 - ™ j+ l/2 + Wj- 1/2 
« = I [ -<- 1/2 -<- 1/2 -<- 1/2 


+<+1/2 + <+1/2 + <+1/2 


= 0 . 


(94) 

(95) 

(96) 


If we define (Af) ■ i — 1,2, to be smoothly positive, then 
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and the additional artificial dissipation operator is given by 


Dl = P-'EpiitAjJlZVf, 

i= 1 

The resulting energy-stable scheme is obtained by adding the additional dis- 
sipation term to the corresponding conventional WENO scheme as 

D 3 = D 3 + D 3 ad . 


8.3 ESWENO Schemes of Seventh- and Eighth-order 


The /■ , l , defined for the seventh- and eighth-order WENO schemes, is 
J'2 



~ W j+l/2fj+l/2 + W j+l/2fj+l/2 + W f+l/2fj+l/2 


PC 


+ w 


R fR 
j+l/2Jj+l/2 


+ W 


RR f HR 
j+l/2J j+1/21 


(97) 


where w LL , w L , w c , w R , and w RR are weight functions that are assigned to 
each respective stencil, and are given by equations (58) and (59) with 


(fj - 3 - 7 fj _ 2 + - 35 fj + 35 f j+1 - 21 f j+2 + 7 f j+3 - f j+4 ) 2 , for <p ± 0 

1 {-fj - 3 + 6/^-2 - 15/j-i + 20 fj - 15/j+i + 6 fj +2 - fj+ 3) 2 , for ip = 0 


The fourth-order linear fluxes /jy , , 2 for r = {LL, L, C, R, RR}, which are used 
in equation (97) are defined as 


f f LL {u j+ 1 / 2 ) ^ 
/ L K+i/2) 
/ C K+i/2) 

/ R (%+i/ 2) 

V / RR K + 1 / 2 ) y 


12 


^ —3 13 —23 25 0 0 0 0 ^ 

0 1 -5 13 3 0 0 0 

0 0 -1 7 7 -1 0 0 

0 0 0 3 13 -5 1 0 

v 0 0 0 0 25 -23 13 -3 } 


f f{uj- 3 ) N 
f(uj- 2 ) 
f(uj- 1 ) 
/(%) 
/K+i) 
f(u j+ 2 ) 
/K+ 3 ) 


(98) 


V 



The derivative matrix D ' is now decomposed into symmetric and skew-symmetric 
parts as 


D 7 


D 


7 

skew 


+ £ 


7 

st/m * 
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As with the fifth-order case, the skew-symmetric component of D 7 takes the 
form 

D 7 s kew = p 1 ( ? 7 ; Qi + Q7 — 0 ; P ~ Axl . 

The scheme converges with seventh-order accuracy if the weights are equal to 
their preferred values 

d LL = P-ip ,d L = g-8 if ,d c = || ,d R = |+8 ip ,d RR = ip (99) 


and with eighth-order accuracy for the specific value ip c = P. 
The D 7 sym can be simplified into the form 


Dl ym = P- 1 (Bi 4 AJ [-Di 4 ] r + -D1 3 A l [£>,+ + D1 2 [B+ 

+ d 1 , a;[d, 1 ] t + B 1 °A 0 7 p 1 °] T ) 


( 100 ) 


where Aj is a diagonal matrix with expressions for the jth element defined by 


(*D* 

1 

" 8 

w j+ 7/2 W f+l/2 

(101) 


1 

24 

[ + W ^5/2 - 9w “ 7 /2 + ^+5/2 

-uf+1/2 + 9^1/2 - ™f+l/2 

(102) 


1 

24 

[ +2 W RR ^i 2 — 11 W RR 5 / 2 + 9w RR 7 / 2 




+ 2w j+3/2 ~ ^ w j+5/2 
- W f+l/2 + ™f+3/2 

(103) 


-!”2 W f-l/2 ^ W j+l/2 

~ 9w f-3/2 + 2 - 2w f+: L/2 


( X l)jJ = YA [ + Gw ftl/2 ~ 13w i+3/2 + 10w i+i 5/2 - 3w i+' 7/2 

+ 3 <+ 1/2 - ^f+3/2 + W i+5/2 
+ W ?-l/2 ~ W f+3/2 
~ W f-3/2 + 4 ^l/2 - 3W&1/2 

+ 3 w i - 1/2 - 10 « h - 1/2 + l ? ylL f R X /2 ~ Gw *«/2 


(104) 
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( 105 ) 



1 

2 


~ W l-\/2 ~ Wi- 1/2 - W ?-l/2 ~ W tl/2 ~ W ?-l/2 

+ W i+ 1/2 + <+ 1/2 + 1/2 + < 5 - 1/2 + W ?+l/2 


0 . 


If we Define (A J) . ., i = 1,4, to be smoothly positive as 

(\ 7 )/,/ = 2 1\/ + ^ _ (« i 
the additional artificial dissipation operator is given by 

D 7 ad = P-'UDAmDif , 

i= 1 

and the resulting energy-stable scheme is obtained by adding the additional 
dissipation term to the corresponding conventional WENO scheme as 

D 7 = D 7 + D 7 ad . 
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